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ABSTRACT 
An  automatic  elevator  controller  for  the  final  phase  of  an 
Instrument  Landing  System  approach  is  designed  using  optimization 
theory  and  the  practicality  of  the  controller  investigated.  The 
problem  is  discussed  and  the  assumptions  stated.  Then  a  mathe- 
matical model  for  the  aircraft  and  a  desired  flare-out  approach 
path  are  derived.  The  aircraft  and  approach  limitations  are  estab- 
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lished  and  the  model  is  teste™.  Dynamic  Tfrogramming  and  the 
Parametric  Expansion  Method  provide  the  optimal  control  from  which 
the  controller  is  designed.  A  computer  program  is  developed  to 
investigate  the  controller.  The  results  are  inconclusive  and  a 
recommendation  for  further  study  is  made. 
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1.   INTRODUCTION 

This  study  is  an  attempt  to  synthesize  an  automatic  control 
mechanism  for  the  terminal  phase  of  an  aircraft  landing  using 
optimization  theory.  The  problem  was  originally  presented  by 
F.  J.  Ellert  and  C.  W.  Merrian  [2],  Dynamic  Programming  [1]  and 
the  Parametric  Expansion  Method  [2]  are  used  to  design  a  simple 
controller. 

Section  two  describes  the  problem  and  presents  the  assumptions 
made  to  simplify  the  initial  study.  Originally  the  investigation  of 
Ellert  and  Merriam  was  to  be  continued  but  after  several  weeks  of 
investigation,  it  was  found  that  their  equations  for  the  aircraft 
represented  only  the  changes  about  an  equilibrium  condition  [3] 
rather  than  the  true  signals.  Equations  representing  only  the  changes 
are  derived  in  Appendix  A  and  are  used  in  section  three  to  obtain 
different  equations  for  the  actual  aircraft  signals. 

In  section  four  a  more  realistic  desired  approach  is  defined 
and  equations  for  the  desired  values  of  the  measured  signals  are 
derived.  The  performance  limits  and  system  constraints  are  dis- 
cussed in  section  five.   Section  six  covers  the  testing  of  the 
aircraft  model. 

Section  seven  develops  the  theory  to  obtain  the  optimal 
control  in  terms  of  the  measured  signals,  aircraft  constants,  and  func- 
tions of  time*  The  necessary  first-order,  linear  differential  equation 
to  obtain  the  time  functions  for  the  optimal  control  are  derived 
in  Appendix  B.   Section  eight  outlines  the  mechanization  of  a 
controller  using  the  optimal  control  equation.   Section  nine 
discusses  the  computer  programming  and  the  investigation 
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conducted  with  the  results  in  section  ten  and  the  conclusions 
in  section  eleven. 
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2.  PROBLEM  STATEMENT 

Landing  an  aircraft  in  inclement  weather  when  the  pilot  has 
no  visual  contact  with  the  runway  until  the  last  few  hundred  feet 
is  a  function  performed  at  many  airports  by  an  Instrument  Landing 
System,  commonly  referred  to  as  ILS. 

The  aircraft  follows  an  electronic  beam  which  is  set  to  a 
glide  path  of  approximately  three  degrees  and  provides  an  azimuth 
reference  to  the  runway.  As  the  aircraft  descends,  the  beam  resol- 
ution becomes  poor  and  sizable  altitude  errors  can  occur.  The  three 
degree  glide  path  is  normally  followed  to  a  point  within  half  a  mile 
from  the  approach  end  of  the  runway.  The  pilot  then  takes  over 
visually  to  land  the  aircraft  or  is  waved  off  to  try  again. 

At  this  time,  a  choice  must  be  made  by  the  pilot  to  follow 
ILS  and  hold  the  rate  of  descent  he  has  or  to  make  a  flare-out 
approach  to  touch  down  on  the  runway  with  a  lesser  rate  of  descent. 
The  flare-out  approach  puts  less  stress  on  the  aircraft  and  is  a 
more  comfortable  landing. 

A  controller  is  needed  to  accomplish  the  flare-out  approach  so 
that  if  visibility  is  restricted,  this  kind  of  landing  may  still  be 
made.  The  following  assumptions  are  made  to  simplify  these  first 
studies: 

a.  The  aircraft  is  laterally  aligned  with  the  runway.  Only 
errors  in  the  vertical  plane  will  be  considered. 

b.  At  one-half  mile  from  the  end  of  the  runway,  the  aircraft 
will  be  waved  off  if  altitude,  rate  of  descent,  pitch  angle,  and 
pitch  rate  exceed  specified  limits. 

c.  Wind  effects  will  not  be  considered. 
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3.   SYSTEM  MODEL  DERIVATION 

The  first  step  in  the  synthesis  of  an  aircraft  control  system 
is  the  development  of  a  mathematical  model  relating  the  control 
(in  this  case,  elevator  deflection)  to  some  measurable  response 
variables.  It  is  essential  that  any  model  chosen  be  tested  to 
insure  that  it  follows  all  known  responses  of  the  real  system.  This 
will  be  discussed  further  in  section  six. 

Figure  1  defines  the  aircraft  coordinates,  angles,  and  positive 
direction  of  elevator  deflection.  Equations  (3.1)  and  (3.2)  have 
been  derived  to  describe  the  motion  of  an  aircraft  [2]  under  the 
following  restrictive  conditions: 

a.  The  equations  are  linearized  with  the  assumption  that 
deviations  from  an  equilibrium  flight  path  are  small.  Only  the 
changes  about  this  equilibrium  condition  are  represented  by 
equations  (3.1)  and  (3.2). 

b„  The  glide  path  angle,  y  ,  is  small  enough  that  the  small 
angle  approximations ,  siny-Y  and  cos  y  =  1  can  be  used.  This 
approximation  is  valid  for  the  landing  paths  to  be  investigated. 

c.  It  is  assumed  that  the  aircraft  velocity,  V,  is  held 
constant. 

Note  that  should  a  controller  be  found  that  is  feasible  from 
this  basic  synthesis 8  a  more  realistic  model  of  the  aircraft  should 
be  used  for  further  studies. 
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Figure  1.  Definition  of  Aircraft  Coordinates,  Angles,  and  Elevator 
Deflection. 

K  W  2(T  8  +  1) 

(3.1)  A9(s)=   S  8   S A6(s) 


s2  +  2rtjs+w 
s  s     s 


(3.2)   Ah(s) 


(T  s  +  1) 
s 


A6 


(s) 


where  A6(s)  =  Laplace  transform  of  the  change  in  pitch  rate,  A6  (t), 

A6  (s)  »  Laplace  transform  of  the  change  in  elevator 

deflection,  A6  (t). 

Ah(s)  ■  Laplace  transform  of  the  change  in  vertical 

acceleration, A  h(t). 

K  i  T  »  w  »  I  a*~e  constants  with  different  values  dependent  on 
s   s   s 

the  type  of  aircraft. 

Equations  (3.1)  and  (3.2)  can  he  combined  to  provide  a  transfer 
function  relating  the  change  in  vertical  acceleration  and  elevator 
deflection, 


15 


2 

(3.3)  Ah(s)  -      s  s  A6^ 

«  +  2rW  s  +  W  2 

s  £       S 

which  can  be  written  as  a  fourth  order  linear  differential  equation 

(3.4)  Ah(t)  +  2^W  Ah(t)  +  W  2Ah(t)   -  K  W  A6(t) 

s        s  s  s 

relating  the  change  in  the  control  variable  to  changes  in  the 
response  variables:   change  in  altitude, A  h(t);  change  in  rate  of 
ascent,  Ah(t);  change  in  vertical  acceleration,  Ah(t);  and  change 
in  jerk,  Ah(t). 

All  of  these  signals  are  not  directly  measurable.   Signals 
that  are  readily  measurable  in  the  aircraft  are  desired.  These  are: 

a.  Altitude,  h(t),  which  can  be  measured  from  a  barometric 
altimeter  or  radar  altimeter. 

b.  Rate  of  ascent,  h(t),  measured  with  a  barometric  rate  meter, 

c.  Pitch  angle,  6(t). 

d.  Pitch  rate,  6(t),  both  measurable  from  gyros. 

These  are  expressed  in  (3.5)  as  an  equilibrium  condition  plus  a 
change  about  that  signal. 

(3.5)  %±(t)  »  hg(t)  +  Ah(t) 
x2(t)  -  h@(t)  +  Ah(t) 
x,(t)  =  e  (t)  +  A9  (t) 
x4(t)  -  ee(t)  +  A6  (t) 

where  the  subscript  e  denotes  the  equilibrium  condition. 

A  set  of  first-order  differential  equations  can  be  developed 
from  (3.2)  and  (3.4)  for  the  changes  about  equilibrium  and  result 
in  equations  of  the  form 

(3.6)  z(t)  »  Az(t)  +  Bu(t) 
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where  js(t)  is  a  4x1  vector  and 

zT(t)  =   [Ah(t),  Ah(t),  A6(t),  A6(t)] 
A  is  a  known  4x4  matrix  of  aircraft  constants 
B  is  a  known  4x1  distribution  vector  of  the  scalar  control, 
u(t). 
The  derivation  of  (3.6)  is  given  in  appendix  A  and  represents 

the  change  of  the  measurable  signals,  x(t),  about  the  equilibrium 

condition.  Equations  (3.5)  can  now  be  written 

(3.7) 


*x(t> 

-    he(t) 

+ 

2x(t) 

x2(t) 

=    he(t) 

+ 

z2(t) 

*3<t) 

■    a.(t) 

+ 

z3(t) 

v> 

-    e.(t) 

+ 

«4(t) 

The  equilibrium  condition,  x  (t),  is  chosen  to  be  a  constant  rate 
of  descent,  as  variations  about  an  ideal  ILS  approach  are  to  be  invest' 
igated.  x  (t)  becomes 

(3.8)  he(t)  -  h(0)  +  h(0)t 

he(t)  =  h(0) 
ee(t)  =  ©(0) 
ee(t)  =  o 

Then  (3.7)  can  be  expressed  as 

(3.9)  xx(t)  =  h(0)  +  h(0)t  +  z2(t)  =  x1(0)  +  x2(0)t  +  Zj(t) 

-  x2(0)  +  z2(t) 


x2(t) 

" 

h(0) 

+ 

z2(t) 

x3(t) 

s 

e(0) 

+ 

z3(t) 

*4<t) 

ss 

*4<t) 

x3(0)  +  z3(t) 


where  z  (0)  =  z2(0)  =  z3(0)  -  0  and  z4(0)  =  x4(0) 
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A  set  of  linear  differential  equations  describing  x(t)  is 

(3.10)  x  (t)  ■  n(0)  +  z  (t) 

X2(t)   =   E2(t) 

x3(t)  -  z3(t) 

X4(t)   =   24(t) 

substituting  (A. 9)  in  (3o10)  gives 

(3.11)  i^t)  =  h(0)  +  z2(t) 

x2(t)  -  a22z2(t)  +  a23z3(t) 

x4(t)  -  a42z2(t)  +  a43z3(t)  +  a44z4(t)  +  b4u(t) 
using  (3 .9) s  (3.11)  Q&n   be  written 

(3.12)  x1(t)  =  x2(t) 

i2(t)  =  a22[x2(t)  -  x2(0)]  +  a23[x3(t)  -  x3(0)] 

x3(t)  -  x4(t) 

x4(t)  -  a42[x2(t)  -  x2(0)]  +  a43[x3(t)  -  x3(0)] 
+  a44x4(t)  +  b4u(t) 
which  are  the  desired  linear  equations  as  a  function  of  the  signals, 
x(t)s  and  the  control,  u(t).  The  system  model  is  also  dependent  on 
the  equilibrium  values  of  the  descent  rate  and  pitch  angle  because 
of  the  chosen  equilibrium  condition. 
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4.   DESIRED  FLARE -OUT  APPROACH 

At  time  t  =  0  with  the  aircraft  some  horizontal  distance,  D, 
from  the  end  of  the  runway,  it  is  desired  to  have  the  aircraft  leave 
the  three  degree  glides  lope,  increase  its  rate  of  descent,  then  flare- 
out  to  arrive  over  the  end  of  the  runway  at  time,  t  =  t. ,  at  a  lower 
altitude  than  if  it  had  remained  on  the  glideslope  beam  and  with  a 
lower  rate  of  descent.  From  time,  t  =  t. ,  until  time  of  touchdown, 
t  =  tf,  the  rate  of  descent  is  held  constant  and  all  accelerations 
are  zero.  Figure  2  shows  this  desired  path.  An  exponential  altitude 
path  over  the  time  interval,  t.  <  t  <  tft   would  be  more  realistic 
but  the  constant  rate  of  descent  was  chosen  for  simplicity. 


ALTITUDE 
(it) 


I 


1(0) « -3 


ILS  GLIDESLOPE 


END  OF 
DISTANCE    (ft)      runway 


t  st, 


TOUCHDOWN 
POINT 


t  *  tf 


trO  u  =  ci 

TIME  Csec) 

Figure  2.  Desired  Flare-out  Approach 
Note:  All  desired  approach  variables  will  be  subscripted,  d,  to 
differentiate  from  the  actual  aircraft  model  signals. 

The  following  assumptions  are  made: 

a.  Origin  of  the  ILS  glideslope  beam  and  the  desired  touchdown 
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point  coincide  1000  feet  down  the  runway. 

b.   Prior  to  t  =  0,  the  aircraft  has  been  following  an  ideal 
ILS  approach. 

From  the  problem  definition  and  the  general  desired  approach 
specifications,  the  known  conditions  must  be  compiled  to  derive 
equations  for  the  desired  measurable  signals. 

First  the  known  conditions  at  times  t  =  0,  t  =  t. ,  and  t  ■  tf 
will  be  stated.  Using  these  boundary  conditions,  equations  will 
then  be  derived  for  the  time  intervals  0  <  t  <  t.  and  t-  <  t  <  tf. 

Prior  to  time  t  ■  0,  the  aircraft  is  descending  on  an  ideal  ILS 
approach.  Since  all  measurable  signals  are  continuous,  the  following 
conditions  exist  at  t  -   0; 

0 

h,(0)9  foj(0),  0,(0),  cu(0),  Yj(0)  are  constants  and 

hd(0)/hd(0),  ed(0),  ad(0),  yd(0)  are  zero. 

The  flight  path  specified  over  the  interval  t.  <  t  <  t-  is 

the  same  as  before  t  ■  0  (constant  rate  of  descent).  Continuity 
requires  the  same  conditions  at  t.  and  tf  with  different  constants. 

The  desired  flare-out  approach  path  over  0  <  t  <  t.  in  terms 
of  altitude,  hrf(t)9  and  derivatives  can  be  approximated  by  a 
seventh  order  polynomial  since  there  are  eight  boundary  conditions, 
four  at  t  =  0  and  four  at  t  =  t. .  These  are 

»  c  c  -f 

(4.1)  hd(t)  -  hd(0)  +  hd(0)t  +  a4t*  +  a^  +  a6t  +  a?t 

3      4      5      6 

h.(t)  ■  hd(0)  +  4a,t  +  5a_t  +  6a,t  +  7a_t 

2       3       4       5 
h*(t)  ■  12a.  t  +  20act  +  30a,t  +  42a_t 
dv         4       5       6       7 

°hd(t)  ■  24a,t  +  60a5t  +  120a,t  +  210a_t 


where  a4  =  [20Y(0)  +  ^(t^)]   V    -  35  fh^O)  -  h^)] 

(tp3    (tp4 

a5  =  [45y(0)  +  397(^)1  -  V    +   84   [hj«»  -  h^)] 

a6  =  [36Y(0)  +  347(^)3   V    -   70   [hd<0)  -  h^)] 

(tp5     (tx)6 

a?  =  [10y(0)  +  KMt^)]  -  V    +   20   [hd(0)  -  h^t^] 

For  the  interval  t.<t<t.,  the  desired  flight  path  is  a 
constant  rate  of  descent.  The  desired  equations  are: 

(4.2)  hd(t)  =  hd(t1)  +  hd(t1)[t-t1] 

hd(t)  =  hd(t]L) 

hd(t)  =  0 

h"d<t)  «  0 

The  desired  pitch  angle,  6  ,(t),  can  be  expressed  as  its  equili- 
brium value,  6d(0),  plus  changes  about 9,(0): 

(4.3)  6d(t)  =   6d(0)  +  A0  d(t). 

Taking  the  derivative  gives 

(4.4)  8d(t)  =   A9d(t). 

From  (3.2), 

(4.5)  A9  (t)   =  _Ts_  Ah"(t)  +   1  A  h  (t). 

V  V 

Integrating  (4.5) 


. . 


(4.6)  A6  (t)   »   Ts  A  h  (t)  +   1  A  h  (t). 
a  y     o        V 

The  desired  values  of  (3.5)  can  now  be  written  in  terms  of  the 
derived  equations 
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(4.7) 


hd(0)  +  hd(0)t  +  a,t4  +  at5  +  a  t**  +  at7, 


x-d(t)  »  <  0  <  t  <  t- 


hjCtj)  +  hd(t1)[t-*1J,.t1  <  t  <  tf 
hd(0)  +  4a,t3  +  5a  t4  +  6a&t5  +  7a_t  , 


x,,d(t)  -  <  0  <  t  <  t 


VV*  t±  <  t  <   tf 


x 


T  0  ^  A  Q 

^ld(0)  +   s  (12a4t  +  20a5t  +  30a6t  +  42a- t) 
(d(t)  ■  ^  +  _!_  (4a,  t  +  5a5t  +  6a&t  +  7a?t  ),  0  <  t  <  ^ 


8d(t1),  tj  <  t  <  tj 


X 


4d 


'  1»  (24a4t  +  60a5t2  +  120a6t3  +  210ayt4) 


,«.s  7+1     (12a. t2  +  20a_t3  +  30a,t4  +  42a,t5) 

(t)     s   \       —  4  5  6  7 


<     t     <     t. 


.  0,    tj     <     t     <     tf 
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5.   PERFORMANCE  LIMITS  AND  SYSTEM  CONSTRAINTS 

Aircraft  constants,  numerical  values  of  the  desired  flight 
path,  and  variable  limits  must  be  specified  before  testing  the 
system  model.  The  T-28  aircraft  was  chosen  for  this  investigation. 
Values  of  the  aircraft  constants  in  (3.1)  and  (3.2)  for  the  T-28 


in 

a 

landing 

configuration 

are: 

(5, 

1) 

K 
s 

■  -1.350  sec"' 

Ts 

■  2.205  sec 

W   ■  0.6161  radians/sec 

s 

g   ■  0.4230 
V  m     157.8  ft. /sec 
The  constants  for  equations  (4.7)  are  calculated  from  (5.1) 
with  reference  to  figure  2.   ILS  weather  minimums  of  visibility 
vary  from  a  half  to  a  quarter  mile.   Using  this  as  a  basis,  a  value 
of  one-half  mile  is  chosen  for  D.   In  consistent  units, 

(5.2)  D  =»  2640  ft. 
Then  from  figure  2, 

(5.3)  hd(0)  m     (D  +  1000)tan  3°  -  190.8  ft. 

Using  the  small  angle  assumption, 

(5.4)  hd(0)  -   yd(0)V 

where  Y.(0)  ■  -0.05236  radians 
and 

(5.5)  t.  -  _D__  sec 

1      V 

tf  -  (D  +  1000)   sec 
V 
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A  value  of  twenty  feet  is  chosen  for  the  desired  altitude 
over  the  end  of  the  runway,  then 

(5.6)  ^d^i)  s  20  ft» 

Yj(t, )  E  =0.02  radians 
d  I 

^d^l*  s  Yd^fcPV  s  -3.156  ft. /sec 

The  physical  limitations  of  the  aircraft  must  be  considered; 
important  details  are  discussed  below. 

1)  Angle  of  Attack.  The  T-28  aircraft  stalls  at  an  angle  of 
attack,  a(t),  of  approximately  21  and  a  speed  of  72  knots.  As 

the  angle  of  attack  approaches  this  limit,  the  speed  cannot  be  main- 
tained constant  and  the  linearized  system  equations  no  longer  repre- 
sent the  aircraft  motion.  It  is  assumed  that  realistic  representation 
is  loet  for  a(t)  >  18°.  This  value  will  be  called  the  stall  angle  of 
attack,  a  .  For  the  T-28  in  landing  configuration  descending  on  an 
ILS  glideslope  at  93.5  knots,  the  angle  of  attack,  a(0),  is  11  or 
0.1920  radians. 

2)  Pitch  Angle.  The  pitch  angle  at  touchdown  must  be  between 
the  limits  0  <  9(t„)  <  14  .  The  lower  limit  prevents  the  nose 
wheel  from  touching  first  and  the  upper  limit  prevents  the  tail  skag 
from  dragging  on  the  runway. 

3)  Elevator  Control.  The  elevator  motion  is  limited  by  the 

mechanical  stops  of  the  actuator  arm.  Assumed  values  are  -35  and 

o 

+15  .  The  direction  is  defined  in  figure  1.  Aircraft  motion  with 

the  elevator  against  these  stops  for  non-zero  time  is  not  permitted 

for  the  linearized  model  used. 
From  figure  1, 

(5.7)  6(t)  ■  a(t)  +  y(t) 
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and  the  desired  pitch  angle  at  t  ■  0  is 

(5.8)  ed(0)  -   ad(0)  +  yd(0)  -  0.1396  radians. 

Since  the  angle  of  attack  must  remain  below  the  stall  value, 
a  ,  an  equation  for  ct(t)  is  desired.  From  (3.2),  the  change  in 

8 

pitch  rate  is 

(5.9)  A9(t)   *   Ts  Ah(t)  +  _1_  Ah(t). 

V  V 

(5.7)  may  be  written  in  the  form 

(5.10)  6  (t)  +   A6(t)   -   a  (t)  +  Aa(t)  +  Y  (t)  +  AY(t). 

e  e  e 

*  •       » 

Taking  the  time  derivative  with  e  (t)  ■  a  (t)  =  Y  (t)  =  °» 

6  6  6 

(5.11)  A9(t)   =   Aa(t)  +  AY(t) 
and  equating  (5.9)  and  (5.11)  gives 

(5.12)  Aa(t)  +  AY(t)   =   Ts   Ah(t)  +  _L_  Ah(t). 

V  V 

It  can  be  shown  [3]  that  Aa(t)  =  T  Ay  so  the  T  term  of  (5.12)  is 

s  s 

associated  with  Aa(t)  and 

(5.13)  Aa(t)  »   T8   Ah"(t) 

V 

(5.14)  AY(t)  =   1  Ah(t). 

V 

Integrating  equations  (5.13)  and  (5.14)  yields 

(5.15)  a(t)   -  a  (t)  +   Ts  Ah(t) 

e         V 

(5.16)  y(t)   =  Y^(t)  +  _l_Ah(t) 

V 

where  a  (t)   =  a(0)  and  Y^(t)   =  y(0). 
e  e 

From  (5.4),  y (0)  ■  n(0)  and  (5.16)  becomes 

V 

(5.17)  Y(t)  =   1  Xg(t) 

V 

Using  (5.7),  a(t)  can  now  be  determined  from  the  system  model 
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by 

(5.18)  a(t)  =  x  (t)  -  _J^_  x  (t) 

In  section  2,  it  was  stated  that  if  the  measurable  signals, 
x(o),  were  outside  specified  limits,  the  aircraft  would  be  waved 
off  for  another  approach.  These  limits  are  somewhat  arbitrary  and 
should  be  chosen  only  after  a  close  study  of  existing  conditions  of 
actual  aircraft  making  ILS  approaches. 

For  this  study  the  following  limits  are  assumed: 

(5.19)  xx(0)  =  xld(°)  ±  20  ft- 


x2(0)  =  x2d(0)  ±  2  ft-/sec 


x3(0)  =  x3d(0)  +  0.0350  radians 


x,(0)  -  x,,(0)  +  0.0043  radians 
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6.  MODEL  TESTING 

Every  mathematical  model  derived  for  a  physical  system  should 
be  tested  by  applying  known  controls  to  insure  that  the  model 
accurately  describes  the  system.  Originally  the  £(t)  equations  were 
used  as  the  aircraft  model  assuming  that  they  represented  the 
signals,  x(t) .  In  this  problem  there  are  two  controls  that  can  be 
used  for  testing.  The  control,  u(t)  =  0,  applied  to  the  model  allows 
the  aircraft  to  follow  its  equilibrium  flight  path  which  is  a  constant 
rate  of  descent  to  the  runway.  From  (3.4)  a  control  can  be  derived 
using  the  equations,  x,(t),  since  they  are  continuous  and  known 
functions  of  time. 

(6.1)  u(t)     =         1  A?(t)     +     2g     Ah .(t)     +       1   A  K.(t) 

wtnT2  KVW  K  V 

K  VW  s     s  s 

s     s 

Since  h  (t)  and  h  (t)  equal  zero, 

(6.2)  /^(t)  =  hd(t) 
A*hd(t)  -  hd(t) 

Using  (4.1)  and  (4.2), 

(6.3)  ^d(fc)  a  l2a4t2  +  20a5t3  +  30a6t4  +  40a?t5 

A*hd(t)  -  24a4t  +  60a5t2  +  120a6t3  +  210a?t4 
A  h^(t)  =  24a4  +  120a5t  +  360a6t2  +  840a?t3 

f or  0  <  t  <  t,. 

(6.4)  Ahd(t)  =  Ahd(t)  -  Ah^(t)  =  0 

for  tx     <     t     <    tf  . 

With  (6.3)  and  (6.4)  used  in  (6.1)  over  their  respective  time  intervals 
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and  this  u(t)  applied  to  the  model,  it  should  follow  the  desired 
flight  path  shown  in  figure  2. 

When  these  controls  were  applied  to  the  first  model,  it  did 
not  follow  either  flight  path  (see  figure  3).  After  deriving  the 
model  in  section  3,  the  same  two  controls  were  applied  and  the 
proper  flight  paths  were  followed  (see  figure  4). 
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7.   APPLICATION  OF  OPTIMIZATION  THEORY 

The  first  step  in  applying  Dynamic  Programming  is  to  obtain 
an  error  index.  This  index  determines  the  form  of  the  resulting 
optimal  control.   The  errors,  x(a)  -  x,(a),  for  0  <  a  <  t-; 
the  errors  at  the  final  time,  x(tf)  -  x,(t_);  and  the  control, 
u(t),  are  considered  to  be  significant  for  this  problem.   The  Bolza 
form  of  error  index  is  considered  adequate  to  insure  satisfactory 
system  performance.  The  error  index,  J,  is  defined  by 

(7.1)  J  =  1  -f[x(tf)  -  xd(tf)]T  F[x(tf)  -  X|<tf)]} 

+  ±J   f  (fi(t)  -  xd(t)]TQ(t)  [x(t)  -  x^t)]  +  u2(t)}  dt 

where  F  is  a  4x4,  positive,  semi-definite,  constant  matrix. 

Q(t)  is  a  4x4,  positive,  semi-definite,  time  variant  matrix. 

Both  F  and  Q(t)  are  arbitrary  weighting  functions  of  the 

system  errors  and  are  chosen  to  be  symmetric. 
The  condition  by  which  (7.1)  is  a  minimum  is  given  by  [1] 

(7.2)  min  \l   f  [x(a)  -  x,(a)]TQ(a)  [x(a)  -  x  (a)]  +  u2(a)}+  dE  j  «  0 
u(a)|_2  l  J   daj 

where  a  is  an  arbitrary  instant  of  time  in  the  interval, 0  <  a  <  tfl 

tf  is  the  fixed  final  time.   E  is  defined  as  the  minimum  error 

function,  E[f(a),  x(a)]  -  min  J,  which  is  assumed  to  be  only  a 

u(a) 

function  of  a  and  the  signals  measured  at  a. 
As  a — ►  tf ,  it  can  be  shown  that 

(7.3)  E[f(tf),  x(tf)]  -  l/[x(tf)  -  xd(tf)]TF[x(tf)  -  xd(tf)]j- 

Define 

(7.4)  H[x(a),  u(a)]  =  l|[x(a)  -  xd(a)]TQ(a) [x(a)  -  xd(*)]  +  u2(a)} 
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then  (7.2)   can  be  written  as 

<7.5)        min     |~H[x(a),   u(a)]  +  djf]  =     0 
u(a)  L  daj 

The  derivative  of  E  with  respect  to  a  may  be  written  as 

(7.6)  dE  =  d  Eff(a),  x(a)1 
da  da 

4 

"  _*E df(a)  +  5T  dxn^Eff(a),  x(a)1 

^f(a)     da      n=l  da      ox  (a) 

The  minimal  error  condition  (7.2)  can  then  be  written 

(7.7)  min  [Hrx(a),  u(a)1  +  hE      df(a)  +  x.(a)pE  + 
u(a)  h  fcf(a)   da  fr 

x2(a)  £E__  +  x3(a)  \E_  +   x4(a)  ^E  1  =  0 
0x2        3x3        ^x4J 

For  (7.7)  to  be  a  minimum  for  all  a,  0  <  a  <  tf,  the  first 
derivative  of  the  term  in  brackets  with  respect  to  u(a)  must  be  zero 
and  the  second  derivative  must  be  positive.  E  was  defined  as  a 
function  of  a  and  x(a)  and  is  not  a  function  of  u(a).  From  the 
system  equations  (3.12)  it  is  seen  that  x,  (a)  is  the  only  time 
derivative  that  is  a  function  of  u(a).  The  derivatives  of  (7.7) 
with  respect  to  u(a)  can  then  be  written 

(7.8)  \ll         +    ftE   )>*4(a)  =  0 
*bu(a)    Tx]T~   ^u(a) 

^2H     +  frE   h2X4(a)    >  0 
bu(a)2     *X4   bu(a)2 

From  (3.12)  and  (7.4) 

(7.9)  b*4(a)   =  b4,   b2^4(a)   =  0 
*U<a>  ou(a)2 

%oH     =  u(a),   fr2H     =  1 
*^  ou(a)2 
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Substitution  of  (7.9)  in  (7.8)  yields 
(7.10)   u(a)  =  -b4  ^E 

1  >  0 

(7.10)  defines  the  optimal  control  for  this  system  over  the  interval, 

0  <  a  <  t  .  Substituting  (7.4),  (3.12),  and  (7.10)  in  (7.7)  gives 

2 

(7.11)  1  j[x(a)  -  xd(a)]T  Q(a)[x(a)  -  x^a)]^  V  [jjE  "    + 

Q.E   df(a)  +  x  (a)  _VE  +  f"a22[x  (a)  -  x2(0)]  + 
"bf(a)  da  bx1    L 

a23[x3(a)  -  x3(0)]l  _^+  x4(a)  JfcE  +  fa42[x2(a)  -  x2(0)] 
J  qx2        bx3   L 

+  a43[*3(a)  "  x3(0)]  +  a44x4(a)]^E  -  ^^M) 

J  ^4       ])*u] 


2 

-  0 


which  is  a  partial  differential  equation  defining  E.   If  (7.11)  can 
be  solved  for  E,  then  the  optimal  control  of  (7.10)  can  be  found. 
The  Parametric  Expansion  Method  [2]  is  now  applied  to  solve 

(7.11)  and  provide  an  optimal  solution.  The  method  is  outlined  as 

follows: 

Assume  a  form  for  E  consistent  with  the  functional  description. 

4  4   4 

(7.12)  E  =  k(a)'  -  2  2  Va)  Xm(a)+2  2f  k  (a)xm(a)x  (a) 

m=l  m=l  p«l  p         p 

where  k  (a)  =  k  (a). 
mp        pm 

Taking  the  derivative  with  respect  to  x4(a),  substitution  in  (7.10) 
yields  the  optimal  control  as  a  function  of  time,  the  aircraft  para- 
meters, and  the  measurable  signals,  x(a). 

(7.13)  u(t)  =  2b4[k4(t)  -  k14(t)xx(t)  -  k24(t)x2(t) 

"  k34(t)x3<t>  *  k44(t)x4(t)1 
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since  a  is  an  arbitrary  instant,  it  can  be  real  time,  t.  Obtain 
the  derivatives  required  by  (7.11)  from  (7.12).   Substitute  these 
in  (7.11)  and  collect  terms  so  that  the  expression  is  in  the  form: 
(7.14)  PQ  +  P1x1(a)  +  P2x2(a)  +  P^a)  +  P^a)  +  P^^a) 
+  P22x2  (a)  +  P33x32(a)  +  P4Ax42(a)  +  P^x^a^a)  + 
P13xx(a)x3(a)  +  P14xx(a)x4(a)  +  P23x2(a)x3(a)  + 

P2AX2(a)x4(a)  +  p34x3(a)x4(a>  "  ° 
The  P  coefficients  contain  terms  of  the  k-parameters,  time 
derivatives  of  the  k(a)'s,  the  desired  signals,  the  weighting  factors, 
the  aircraft  constants  and  the  chosen  equilibrium  conditions.  For 

(7.14)  to  be  valid  for  all  values  of  the  measured  signals,  each  P 
coefficient  must  independently  equal  zero. 

(7.15)  P   -  P   -  P    =0 

o      m      mp 

where  m  ■  1,  2,  3,  4;  p  ■  1,  2,  3,  4. 

Equation  (7.15)  gives  a  set  of  15  independent,  ordinary,  first- 
order  differential  equations  which  define  the  k(a)'s.  These  equations 
with  the  boundary  conditions  are  derived  in  Appendix  B. 


34 


8.   IMPLEMENTATION  OF  A  CONTROLLER 

Several  types  of  controllers  could  be  designed  depending  on 
the  assumed  form  of  E  and  the  methods  available  for  measuring  time 
and  velocity.  The  block  diagram  shown  in  figure  5  is  obtained  from 
(7.13)  and  is  a  simple  analog  control  system.  A  constant  speed  motor 
rotates  the  shaft  at  a  constant  velocity.  The  k(t)'s  of  (7.13) 
represent  the  settings  of  wire-wound  potentiometers  where  the  wiper 
arms  rotate  with  the  shaft.  The  measured  signals  are  multiplied 
by  the  appropriate  k(t),  summed  and  multiplied  to  provide  the  ele- 
vator deflection  signal. 


V*i 


jytjfiw 


Vow) 


■Wflft) 


EACH  POT  IS  WIRE -WOUND 
AS  A  FUNCTION  OF  THE 
RESPECTIVE  k(t). 


Figure  5.  Elevator  Controller 

In  the  derivation  of  the  k(t)  equations,  the  equilibrium  values 
of  the  descent  rate,  x.  (t),  and  pitch  angle,  x«je(t),  must  be  known. 
This  means  that  the  optimal  k(t)  equations  are  dependent  on  the 
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equilibrium  conditions.  This  is  a  function  of  the  aircraft  model 
and  would  be  altered  with  a  different  aircraft  model.  Little  error 
should  be  introduced  if  these  values  are  assumed  to  be  the  desired 
descent  rate  and  pitch  angle  at  problem  start. 
(8.1)   x2e(t)  =  x2d(0) 

x3e(t)  =  X3d(0> 
By  assuming  these  values  the  system  is  no  longer  optimal,  but  it  is 
desired  to  sacrifice  optimality  in  the  interest  of  simplicity  and  to 
investigate  the  feasibility  of  this  simple  controller.  The  k(t) 
equations  are  also  dependent  on  the  aircraft  constants  so  a  diff- 
erent set  must  be  obtained  for  each  type  of  aircraft. 

At  this  point  two  big  questions  remain.   1)  Can  k,(t),  k.,(t), 
k^,(t),  k~,(t),  and  k»,(t)  be  found  such  that,  for  all  aircraft 
with  x(0)  within  the  signal  limits  given  by  (5.19),  the  aircraft 
will  make  a  satisfactory  landing?  2)  Are  the  k(t)'s  found  above 
linear  enough  to  be  wound  on  potentiometers?  A  method  to  investi- 
gate the  answers  to  these  questions  is  outlined  in  section  9. 
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9.   INVESTIGATION  PROCEDURE 

Since  all  of  the  constants  and  variables  of  the  fc(t)  equations, 
(B.5),  have  been  calculated  or  are  known  except  the  weighting  factors 
F  and  Q(t),  the  answers  to  the  questions  in  section  8  can  be  investi- 
gated by  the  following  procedure: 

1)  Arbitrarily  pick  the  weighting  factors  of  the  error  index, 
F  and  Q(t). 

2)  Integrate  equations  (B.5)  from  t  a  t,  to  t  ■  0  to  obtain 
the  k(t)'s  of  (7.13). 

3)  Integrate  the  aircraft  model  equations  (3.12)  from 

t  ■  0  to  t  ■  t-  using  u(t)  of  (7.13)  with  the  k(t)'s  obtained  in  2. 

4)  Observe  the  final  signal  errors,  x(tf)  -  x  (tf).   If  these 
are  not  within  prescribed  limits,  change  F  and/or  Q(t)  and  repeat 
steps  2  through  4. 

A  program  was  written  to  accomplish  the  above  procedure  using  the 
CDC  1604  computer  and  an  existing  modified  Runge-Kutta  integration 
subroutine.  There  were  two  main  sections  in  the  computer  program. 
1)  Equations  (B.5)  were  integrated  from  t=t.  to  t  =  0  with  the 
values  of  k4(t),  k14(t),  k£4(t),  ^(t),  and  k^(t)  stored  in  a 
buffer.   2)  Three  sets  of  equations  (3.12)  were  integrated  from 
t  -  0  to  t  -  t£  using  u(t)  of  (7.13)  and  the  k(t)'s  that  were 
stored  in  the  buffer. 

Three  aircraft  models  were  chosen  for  each  program  so  that  one 
could  be  started  with  positive  errors,  one  started  on  the  desired 
trajectory,  and  the  third  started  with  negative  errors.  This  allowed 
all  of  the  system  errors  to  be  observed  in  one  program. 
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The  final  error  limits  were  arbitrarily  chosen  as: 

a)  An  altitude,  descent  rate  combination  such  that  all  three 
aircraft  land  within  +500  feet  of  the  desired  touchdown  point  with 
a  descent  rate  of  +  100  ft./min  or  +  1.66  ft. /sec  from  the  desired 
rate  of  descent,  XjjCtf)* 

b)  Pitch  angle  as  stated  in  section  5,  0  <  0(tf)  <  14  . 

c)  No  limit  was  set  on  the  pitch  rate  but  it  must  be  small 
for  the  approach  to  closely  follow  the  desired  flight  path. 

For  this  investigation,  the  aircraft  was  considered  to  have 
only  an  altitude  error  and  Q(t)  ■  Q,  a  constant  matrix  consisting 
of  the  diagonal  elements  only.  The  study  was  conducted  in  the 
following  manner: 

1)  Determine  if  successful  approaches  can  be  made  with  con- 
stant weighting  of  any  one  signal;  or  any  two  signals;  or  combin- 
ations of  any  three  signals;  and  finally,  combinations  of  all 
four  signals. 

2)  Determine  if  successful  approaches  can  be  made  with  constant 
weighting  of  all  signals,  progressively  adding  diagonal  elements  of 
the  F  matrix. 

3)  When  successful  approaches  are  found,  determine  the  nature 
of  the  k(t)'s  required  for  the  controller. 
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10.   RESULTS 

Digital  simulation  was  used  to  study  the  areas  outlined  in 
section  9.  The  results  were: 

1)  Constant  weighting  of  all  four  signals  was  required  to 
obtain  a  satisfactory  trajectory.  All  final  error  limits  were 
met  except  that  of  the  minimum  descent  rate.  Figure  6  shows  the 
altitude  vs.  time/distance  for  three  aircraft  started  with  an 
initial  altitude  error  of  +  10  feet,  where  q...   =  2.5xlO~   , 

q  2  =  8.5x10"   ,  q33  -  6.0xl0*5,  and  q,,   -  8.8x10  .  All 
other  q's  and  F  were  zero.  The  minimum  descent  rate  at  touchdown 
was  -5.36  ft. /sec  compared  with  the  limit  of  -4.91  ft. /sec. 

2)  This  area  was  not  fully  investigated  but  successful 
approaches  were  made  with  constant  weighting  of  the  four  signals 
(Q  equal  to  zero  except  q...,  q22>  q~,,  q,,)  and  weighting  all 
four  signals  at  the  final  time  (F  equal  to  zero  except  f--,  f~~, 
f~,  £/./,)•     Figure  7  shows  the  altitude  vs.  time/distance  for  three 
aircraft  started  with  an  initial  altitude  error  of  +  20  feet  where 
qn  -  7.10X10"1?  q22  =  2.95xl0~7,  q33  -  7.00xl0*5, 

q44  -  8.00xl0"4  and  fn  -   1.80xl0~9,  f22  -  2.00xl0~7, 

-3  -3 

f00  -  2.90x10  ,  £,,      -  3.00x10  .  All  final  errors  were  within 
33  '  44 

specified  limits. 

3)  For  the  approach  of  part  1,  where  F  is  zero,  k, (t)  is  shown 
in  figure  8,  k34(t)  in  figure  9,  and  ^(t)  in  figure  10.  ^(t)  and 
k«,(t)  were  insignificant- (i.e.  nearly  zero). 

For  the  approach  of  part  2,  k,(t)  is  shown  in  figure  11,  ^.(t) 
in  figure  12,  and  k^Ct)  in  figure  13.  Again  ^14(t)  and  k2A(t)  were 
insignificant. 
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11.   CONCLUSIONS 

For  the  results  specified  in  section  10,  the  following  conclu- 
sions are  stated: 

1)  Although  the  minimum  descent  rate  error  was  not  within 
limits  using  the  Q  matrix  only,  it  is  felt  that  with  several  small 
adjustments  of  the  weighting  factors  the  descent  rate  could  be  brought 
within  limits.  The  k(t)'s  are  of  the  same  sign  with  maximum  values 
less  than  one  and  could  be  wound  on  the  potentiometers. 

2)  While  all  of  the  error  limits  were  met  for  the  approach 
using  both  F  and  Q,  the  k(t)'s  changed  sign  and  are  unsatisfactory 
for  this  controller.  As  adjustments  were  made  to  F  and  Q,  bringing 
the  flight  paths  closer  to  the  desired  trajectory,  the  k(t)'s 
approached  the  same  sign  over  the  interval  of  flight.  It  is  believed 
that  if  F  and  Q  are  adjusted  for  flight  paths  with  smaller  errors, 
this  problem  would  be  eliminated. 

Two  recommendations  are  made: 

1)  That  investigation  be  continued  until  definite  conclusions 
can  be  obtained  about  the  practicality  of  this  controller. 

2)  That  a  method  be  investigated  to  systemmatically  adjust 
the  weighting  factors,  the  elements  of  the  Q  and  F  matrices. 


48 


BIBLIOGRAPHY 

1.  R.  E.  Bellman.   Dynamic  Programming,   Princeton  University  Press, 

Princeton,  New  Jersey,  1957. 

2.  F.  J.  Ellert  and  C.  W.  Merriam.   "Synthesis  of  Feedback  Controls 

Using  Optimization  Theory  -  An  Example,"  IEEE  Transactions 
on  Automatic  Control,  v.  Ill,  April  1963:  89-103. 

3.  CD.  Perkins  and  R.  E.  Hage.   Airplane  Performance  Stability 

and  Control.   John  Wiley  and  Sons,  Inc.,  New  York,  1949. 


49 


APPENDIX  A.  DERIVATION  OF  s(t)  EQUATIONS 
Define 


(A.l)   Z]L(t)  =  Ah(t] 


z2(t)  -  Ah(t) 

z3(t)  ■  A9(t) 

z4(t)  -  &4(t) 

u(t)  -  AS(t) 

The  time  derivatives  of 

the  z(t)'s  are: 

(A. 2)   zx(t)  -  An(t) 

■  «2<tj 

z2(t)  -  dh(t) 

z3(t)  -  A9(t) 

-  «4<o 

z,(t)  -  Ae(t) 

••      • 

Equations  for  Ah(t)  and  A9(t)  in  terms  of  the  states,  _z(t),  are 
needed.   (3.2)  can  be  written  as 


(A. 3)   T4h(s)  +  Ah(s)  =  VA9<s) 
s 

Taking  the  inverse  Laplace  transform  and  solving  for  Ah(t)  yields 

(A. 4)   A*h(t)  m     -_l^h(t)  +  _VA9(t) 

T  T 

s  s 

Differentiating  (A. 4)  gives 

W  •••  • 

(A. 5)      Ah   (t)     ^     -JAh(t)     +    _V,A9(t) 

T  T 

s  s 

•  •  • 

Integrating   (A. 4)   and  assuming  Ah(0)     =     0  gives 

(A. 6)      Ah\t)     -     -_ldh(t)     +    _VA9(t) 

s  s 

Substituting  equations  (A. 4),  (A. 5)  and  (A. 6)  in  (3.2)  and  solving 
for  A9(t)  gives 


(A. 7) 

A9(t)      =     1 
V 

1     -   2SWS 
m  2         T 
J. 

-  r  i   -  2^  +  w 2 

+     Ws^(t) 

+  Ws2lAh(t) 


-   2SW 


sl  AO(t) 


Inserting  (A. 6)  and  (A. 7),  written  in  terms  of  (A.l),  into  (A. 2) 

yields 

(A. 8)   ix(t)  =  z2(t) 

z9(t)  -  -  JL  z9(t)  +  _V  z.(t) 


*3<t>     ■ 

■     «4(t) 

z4(t)     -     1    f~   1     -   2$Wg  +  Wr2 

z2(t)   - 

r                   3 

1        -   2$Wfi     +  W 

«3<t) 

v  LTs     Ts 

T  2 
s             s 

+  111  "  2?wl    z4(t)     +     KsT8Wfi2u(t) 

which  is  in 

the  form 

(A. 9)        z(t: 

1     -     Az(t)  +Bu(t) 

where 

0 

10       0 

"o 

A     s 

0 

a22  a23  ° 

B     s 
1 

0 

0 

0       0       1 

0 

and  a22 

0 

a42  a43  a44 

_     1         q               a 

=r'    23 

T 

s 

V'  ***    "    5 

T                        V 

s 

A 

l 

T2 
s 

-  2$Ws     +     W 

T 

s                   J 

a4 

3     = 

- 

1 

m    2 

T 
s 

2$WS  +  We 
s 

Ts 

,   a 

44     a 

1 

1    -   2$W 

r 

s 

t 

b.=K  T 
4     s 

s  s 
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APPENDIX  B,  DERIVATION  OF  k   (t)  EQUATIONS 

rap 

The  partial  derivatives  of  E  defined  in  (7.12)  are  needed  for 
substitution  in  (7.11).   Since  a  is  an  arbitrary  instant  of  time  in 
the  interval  0  <  a  <  tf,  it  can  be  chosen  as  real  time,  t.  The 
functional  notation^  (t),  on  all  the  variables  will  be  dropped  for 
simplicity. 

The  required  partial  derivatives  of  E  are: 


(B»1)   A3L_  df(t)  "  H  i  +  2^  k   +  dE  k 
U(t)   dt      "fck       5fcT       *k2 

+  IL.  *3  +  M_  \  +  ^S_  kn  +  *E__  Ic22  +  &E__  k33 
ik3  tfk,  dk^  bk22  "&k33 

+  *£_  k«  +  JJ_  k12  +  IS-  kl3  +  Jl£_  ku  +  Jkg_   k23 

^44  »ki2  ^k13  >ku  *k23 

+    >E     k^,  +     >E     k0/ 

Sk         24       Xk         34' 
0R24  *   34 

|E     =      1»  il    "      -2xt»  il     =     "2x9»  il     "      -2xv 

ik  >k1  l     ^k2  *     ^ 

2  2  2 

\E     -     -2x49  \E     -     xx   ,  >>E     ■     x2   ,  ^JL     =     x3   , 

>k4  un  ^k22  W^3 

ay 

\E       -     x.    ,  ^E       -     2x^2,  ^E_     ■     2xxx3,  j^E_     =     2xA» 


*k44  ^2                       *k13                       *k14 

^E  -     2x^x„,  tE       ■     2x-x.,  *E       -     2x-x. , 

*k23  6*24                       *   34 

)>E  -     2(-k1  +  ^j^  +  k12x2  +  k13x3  +  k14*4), 
1 

|E_  »     2C-k2+k12x1+k22x2+k23x3+k24x4), 
**2 


qE_       ■     2(-k,  +  k^x,  +  k0,x,,  +  k^x,,  +  k,AxA), 
ix3 


l3         13  1         23  2         33  3         34  4' 
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IJL    "     2(-k4  +  kux1+  k24x2  +  k34x3  +  k44x4), 
4 


12 

4(k/   -  2k.k_.x_    -  2k.k0/x0   -2k. k«. 


LxJ       "     4(k4 


,.    ,  J.t-22_,i-22_,i2     2J_122 

~  ZK4K4AX4  +  ki4  xi    +  K24  x2    +  k34  x3  +  k44  X4 

+  2ki4k24XlX2  +  2k14k34XlX3  +  2k14k44XlX4 

+  2k24k34X2X3  +  2k24k44X2X4  +  2k34k44X3X4^ 
The  derivatives  of  (B.l)   are  then  substituted  in  (7.11)   to  obtain 


(B.2)     1   f(x  -  xd)TQ(x  -  xd)j 


-2b42(k42  -  2k4k14Xl  -2k4k24x2 


"2k4k34x3  "  2kkk44x4  +  k14  V  +  ht>W  +  Su^l  + 

2  2 
k44  x4  +  2k14k24xlx2  +  2k14k34xlx3  +  2k14k44xlx4  + 

2k24k34X2X3  +  2k24k44X2x4  +  2k34k44X3X4)  +  k  "  2kixi  " 

2k2x2  -  2k3x3  -  2k4x4  +  kuXl2  +  k22x22  +  k33x32  + 

•    2    •         •  •  * 

k44X4  +  2k12XlX2  +  2k13XlX3  +  2k14XlX4  +  2k23X2X3  + 

2k24X2X4  +  2k34X3X4  +  2x2('kl  +  kllXl  +  k12X2  +  k13X3  + 
kux4)  +  2[a22(x2  -  x2(0))  +  a23(x3  -  x3(0))]  (-k2  + 

k12Xl  +  k22X2  +  k23X3  +  k24X4>  +  2x4("k3  +  k13Xl  + 
k23x2  +  k33x3  +  k34x4)  +  2[a42(x2  -  x2(0))  + 

a43(x3  '  X3(0))  +  a44X4]  (~k4  +  k14Xl  +  k24X2  + 

^a  +  W  "  ° 

The  matrices  in  the  first  term  are  expanded  and  all  terms  are 
grouped  in  the  form  of  (7.14).  The  resulting  P  coefficients  are: 
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(B.4) 


P   s  k  - 
o 


2  2 
k/  +  11 


I2  +  ^22X2d2  +  <*33X3d2  + 


q44x4d  }  +  q12Xldx2d  +  q13xldx3d  +  q14xldx4d  +  q23x2dx3d  + 
q24X2dx4d  +  q34X3dx4d  +  2k2(a22x2(0)  +  a23x3(0))  +  2k4(a42x2(0) 


+  a43x3 


(O)) 


Pl  s  "2kl  +  4b4\k14  ■  (qllXld  +  q12x2d  +  q13x3d  + 
q14x4d>  "  2kl2(a22X2(0)  +  a23x3(0))   =  2k14(a42x2(0)  + 


a43x3 


(0)) 


P2  -   »2k2  +  4b4  k4k24  -  (q12xld  +  q22x2d  +  q23X3d  + 
x4d)  -  2k22(a22x2(0)  +  a23x3(0))   -  2k24(a42x2<0)  + 


a43X3 


(0)) 


2kl  "  2a22k2  "  2a42k4 


2 

2k3  +  4te4  k4k34  "  (q13Xld  +  q23X2d  +  q33X3d  + 

"  2k23(a22X2(0)  +  a23x3(0^  "  2k34(a42X2(0)  + 

a43x3(07     ~  2a23k2  "  2a43k4 

P4     -      -2k4  +  ^b42k4k^   -    (q14xld  +  q24X2d  +  q34X3d  + 

<Ux4d>   =  2k24(a22x2(0)  +  a23x3(0>)  "  2k44(a42x2(0)  + 


a43x3 


'3  "  2a44k4 


11     S     kll  +  |qll  "  2\\h 


K^  +  I<Ir 


33 


l\\k    +  2k12  +  2a22k22  +  2a42k24 


2       2 
'4  k34     +  2a23k23  +  2a43k34 


P    ■  k   + 
r44     44 


"  2b4  k44  +  2k34  +  2a44k44 
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P12  »  2k12  +  q12  -  4^^  +  2ku  +  2a22k12  +  2a42ku 

P13  -  2k13  +  q13  "  4b42k14k34  +  2a23k12  +  2a43k14 

P14  "  2k14  +  q14  "  4b42k14k44  +  2k13  +  2a44k14 

P23  B  2k23  +  q23  *  4b4%4k34  +  2k13  +  2a22k23  +  2a23k22 

+  2a42k34  +  2a43k24 
P24  "  2k24  +  ^24  *  4b42k24k44  +  2k14  +  2k23  +  2a22k24 

+  2a42k44+  2a44k24 

P34  "  2k34  +  q34  ■  4b42k34k44  +  2k33  +  2a23k24+  2a43k44 
+  2a44k34 

For  each  P  coefficient  to  equal  zero  for  all  values  of  the 
measured  signals  as  required  by  (7.15),  the  following  first  order, 
linear  differential  equations  must  be  satisfied: 
(B.5) 

k  '  2b4V  "  -l^ll^d2  +  «22x2d2  +  «33x3d2+  ^44x4d2) 

"  q12Xldx2d  "  q13Xldx3d  "  q14Xldx4d  "  q23x2dx3d 

"  q24x2dx4d  "  q34x3dx4d  "  2k2(a223C2(0)  +  a23x3(0)) 

"  2k4(a42x2(0)  +  a43x3(0)J 
kl  =  2\\klU   "  -|(qllxld  +  q12x2d  +  q13x3d  +  q14x4d> 

*  kl2(a22x2(0)  +  a23x3(0>)  "  k14(a42X2<0)  +  a43x3<°>) 
k2  *  2b42k4k24  "  ^(q12xld  +  q22x2d  +  q23x3d  +  q24x4d> 

-k22(a22x2(0)  +  a23x3(0))   "  k24(a42x2(0)  +  a43x3(0)) 

"  kl  "  a22k2  "  a42k4 
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k3     -      2b4\k34   -  I(qoXld  +  q23X2d  +  q33X3d  +   q     x^) 

2 

"  k23(a22x2(0)  +  a23x3(0))     "  k34(a42x2(0)  +  a43x3(0)) 

-  a23k2     -  a43k4 

K    m     2b4\k44  "  |(qX4xld  +  q24x2d  +  q34x3d  +  q44x4d> 

"  k24(a22x2(0)  +  a23x3(0))      "  k44(a42x2(0)  +  a43x3(0)) 

"  k3  "  a44k4 

kll     e     2\\b     ■  |  *11 

z 

k22  2b4^24     "  |  q22  "  2k12  "  2a22k22  "  2a42k24 

'  2   2 

K33       4  34    ^  q33     23  23     43  34 

k44  a  2b42k442  °  |  q44  "  2k34  "  2a44k44 

•  2. 

k12     2b4  14k24  "  |  q12  "  kll  "  a22k12  "  a42k14 

k13  "  2b4  14k34  "  |  q13  "  tt23k12  "  a43k14 

2 

•  2 

k14  B  2b4  k14k44  "  |  q14  "  k13  "  a44k14 

2 

k23  s  ^4^24^4  "  |  q23  *  k13  "  a22k23  "  a23k22*  a42k34 

"  a43k24 
k24     2b4  k24k44  '  |  q24  "  k14  "  k23  "a22k24  "  a42k44  "  a44k24 

k34  s  2b4^34k44  "  |  q34  "  k33  "  a23k24  "  a43k44  "  a44k34 

2 

where  the  boundary  conditions  at  t  s  tf  must  be  specified. 

Equation  (7.3)  is  the  expression  for  E  at  t  ■  t-.  Equating 
this  to  (7.12)  evaluated  at  t  ■  t-  gives 
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4  4         4 

(B.6)       k(t  )   -  2%   k  (t  )  x  (t  )     +5]       gk     (t  )x(t  )x  (t  ) 

t  ,      m     i       m     i  ^^.  ,   mp     r     m     r     p     r 

I  {  te<tf)   -  2Ed(tf)]T  F   [x(tf)   -  Xd(tf)]j 

Expanding  both  sides  and  solving  for  the  k(t.)'s  gives  the  boundary 
conditions  on  the  k's: 

(B.7)        k<t£)     -     I   [£uXld2(t£)  +  f22KM2(tf)  +  f33x3d2(t£)  + 

Vw^V  +  £12xld<t£>x2d(t£>  +  £13xld(tf)x3d(tf)+ 
f14xld<t£>x4d<t£>  +  f23x2d(tf)x3d(t£)  + 
f24x2d<tf)x4d(t£)  +  £34x3d(tf)x4d(,:£) 
W     =    \  [£uKld(t£)  +  f12xM(t£)  +  fuXjjOjP  +  f14"4d(t£)] 

k2(t£>  "  \  [£12xld(t£)  +  f22x2d(t£)  +  f23x3d(t£)  +  f24«M(tf)] 

W  *  |  [£13xld(t£)  +  f23x2d(t£)  +  £33x3d(tf)  +  ^^(tj)] 

W  =  \  [fuxld(t£)  +  £24x2d(t£)  +  f34x3d(t£)  +  V4d(tf)] 

kU<V  "  ifll 

k22(tf>  =  \  f22 

k33(tP  '  |  £33 

k44<tJE)  -  |f44 

k12<*P  =  |f12 

k13<tf)  "  -1  £X3 

k14<Cf>  "  |f14 
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k23<*f> 


k24(tf>     *    |f24 
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